Data from 1960 to 2010 raining days, there are 4168 observations 13 days of which are rain exceeds 1000 (10cm)
Change of exceeds 1000 is 0.003119002 0.31% change
install.packages("moments") # Install moments package
trying URL 'https://cran.ma.imperial.ac.uk/bin/macosx/contrib/4.0/moments_0.14.tgz'
Content type 'application/x-gzip' length 54265 bytes (52 KB)
==================================================
downloaded 52 KB
The downloaded binary packages are in
/var/folders/4x/m08khw3x0497mht_rr4wwxmh0000gn/T//Rtmpmqs2am/downloaded_packages
Data
read.csv("./raw/prcp.csv") %>%
mutate(date = as.Date(date)) %>%
mutate(year = format(date, "%Y")) %>%
filter(prcp != 0) %>%
dplyr::select(date, year, prcp) -> data
data %>%
filter(year <= 2010) %>%
pull(prcp) %>%
hist(breaks = 1000, plot = FALSE) -> hist
Modling Distribution
Use data from 1960 - 2010 for modeling precipitation
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
model <- drm(data = dens, y ~ log(x), fct = EXD.2())
summary(model)
Model fitted: Exponential decay with lower limit at 0 (2 parms)
Parameter estimates:
Estimate Std. Error t-value p-value
d:(Intercept) 0.15567499 0.00099172 156.97 < 2.2e-16 ***
e:(Intercept) 1.27424832 0.00578594 220.23 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error:
0.001088457 (945 degrees of freedom)
d = 0.15567499
e = 1.27424832
k <- 1/e/d
plot(x = dens$x, y = fitted(model), type = "l", col = "blue", xlim = c(0, 300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))
mod.exd <- drm(y ~ x, fct = EXD.2(), data = dens)
plot(x = dens$x, y = fitted(mod.exd), type = "l", col = "blue", xlim = c(0,300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))
summary(mod.exd)
1 - pexp(log(1000)*k, d)
1/d
data %>%
filter(year <= 2010) %>%
pull(prcp) %>%
mean() %>%
log()*k
0.31% change versus 0.442%
Labs:
A Single Years
yr = 1960
data %>%
filter((year >= yr) & (year <= yr + 9)) %>%
pull(prcp) %>%
hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)
yr = 1970
data %>%
filter((year >= yr) & (year <= yr + 9)) %>%
pull(prcp) %>%
hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)
yr = 1980
data %>%
filter((year >= yr) & (year <= yr + 9)) %>%
pull(prcp) %>%
hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$density
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)
Risk Trends
Risks Through out Years
mtr <- data.frame(year = YEAR)
list <- c()
for (i in YEAR){
data %>%
filter((year >= i) & (year <= i + 9)) %>%
pull(prcp) %>%
hist(breaks = 400, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
p <- 1 - pexp(log(1000)*k, d)
list <- c(list, p)
}
plot(1:length(list), list)
Use smaller breaks
Find the right break numbers
YEAR <- c(1960, 1980, 2000)
BREAKS <- c(19, 20, 21, 22)
lenlist <- c()
kulist <- c()
for(yr in YEAR){
data %>%
filter((year >= yr) & (year <= yr + 19)) %>%
dplyr::select(year, prcp) %>%
pull(prcp) %>%
length() -> n
lenlist <- c(lenlist, n)
data %>%
filter((year >= yr) & (year <= yr + 19)) %>%
dplyr::select(year, prcp) %>%
pull(prcp) %>%
kurtosis() -> ku
kulist <- c(kulist, ku)
}
Square_root = c()
Sturges = c()
Rice = c()
Donna = c()
for(n in lenlist){
sq = ceiling(n^0.5)
sg = ceiling(log(n, 2)) + 1
rc = ceiling(2*n^(1/3))
igh = (6*(n - 2)/(n +1)/(n +3))^0.5
dn = log(abs(1 + kulist[which(lenlist %in% n)]/igh) , 2) + log(n, 2) + 1
Square_root = c(Square_root, sq)
Sturges = c(Sturges, sg)
Rice = c(Rice, rc)
Donna = c(Donna, dn)
}
data.frame( Year = YEAR,
Length = lenlist,
Square_root = Square_root,
Sturges = Sturges,
Rice = Rice,
Donna = Donna
) -> meth
meth
Work Flow
In a 10 years interval compute risks of precipitation higher than 10cm (prcp > 1000) by using different breaks of histogram from usgin 250 to 1000
Choosing bin have many choices actually
YEAR <- c(1960, 1980, 2000)
BREAKS <- c(19, 20, 21, 22, 30, 40, 50, 60, 70, 80, 90)
for(n in BREAKS)
data %>%
filter((year >= 1960) & (year <= 1969)) %>%
pull(prcp) %>%
hist(breaks = n)
kurtosis(data$prcp)
list <- c()
listB <- c()
Year <- c()
dlist <- c()
elist <- c()
klist <- c()
for(n in BREAKS){
for(i in YEAR){
data %>%
filter((year >= i) & (year <= i + 19)) %>%
pull(prcp) %>%
hist(breaks = n, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$density
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
p <- 1 - pexp(log(1500)*k, d)
list <- c(list, p)
listB <- c(listB, n)
Year <- c(Year, i)
dlist <- c(dlist, d)
elist <- c(elist, e)
klist <- c(klist, k)
}}
mtr <- data.frame(year = Year,
Breaks = listB,
Probability = list,
parameter_d = dlist,
parameter_e = elist,
parameter_k = klist)
Visulisation
library(wesanderson)
mtr %>%
mutate(Breaks = as.factor(Breaks)) %>%
ggplot(aes(x = year, y = Probability, color = Breaks)) +
geom_point() +
geom_line() +
scale_color_brewer(palette = "Blues") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
) +
ylab("Risks of Rain Exceeds 10cm") +
ggtitle("Risks of Rian Exceeding 15 cm increase evey 2 decades") +
xlab("Two Decades Time From _ ") -> g
#library(plotly)
ggplotly(g)
n too large, allowed maximum for palette Blues is 9
Returning the palette you asked for with that many colors
mtr %>%
mutate(year = as.factor(year)) %>%
mutate(Breaks = as.factor(Breaks)) %>%
ggplot(aes(y = parameter_d, x = parameter_e, color = year)) +
geom_jitter(alpha = 0.8) +
scale_color_brewer(palette = "Reds") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
) +
xlab("parameter e: skewness of the curve") +
ylab("parameter d: height of the curve")

NA
Parameter d influence the height of pdf, parameter e influence the lenghth of pdf… this suggest in general, precipitation distribution cure are getting flatter and
plot_ly(x= mtr$parameter_e, y= mtr$parameter_d, z=mtr$Breaks, type="scatter3d", mode="markers", color=mtr$year)
mtr %>%
dplyr::select(year, Breaks, parameter_d, parameter_e) %>%
mutate(Breaks = as.factor(Breaks)) %>%
melt(id = c("year", "Breaks"), variable.name = "parameter") %>%
ggplot(aes(x = year, y = value)) +
geom_line(aes(x = year, y = value, linetype = parameter, color = Breaks)) +
scale_color_brewer(palette = "Purples") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
)

mtr %>%
mutate(Breaks = as.factor(Breaks)) %>%
ggplot(aes(x = year, y = parameter_k, color = Breaks)) +
geom_point() +
geom_line() +
scale_color_brewer(palette = "Blues") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
)
calculate different
list <- c()
listB <- c()
Year <- c()
dlist <- c()
elist <- c()
for(n in BREAKS){
for(i in YEAR){
data %>%
filter((year >= i) & (year <= i + 9)) %>%
pull(prcp) %>%
hist(breaks = n, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
p <- 1 - pexp(log(600)*k, d)
list <- c(list, p)
listB <- c(listB, n)
Year <- c(Year, i)
dlist <- c(dlist, d)
elist <- c(elist, e)
}}
mtr <- data.frame(year = Year,
Breaks = listB,
Probability = list,
parameter_d = dlist,
parameter_e = elist)
mtr %>%
mutate(Breaks = as.factor(Breaks)) %>%
ggplot(aes(x = year, y = Probability, color = Breaks)) +
geom_point() +
geom_line() +
scale_color_brewer(palette = "Blues") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
)

g <- df %>%
ggplot(aes(x = date, y = prcp)) +
geom_line(color = "blue") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
) +
geom_line(aes(x = date, y = rep(1500, length(date))), color = "lightblue") +
ylab("Precipitation (in tenth of mm")
ggplotly(g)
data %>%
filter(prcp >= 1000)
mtr %>%
filter(Breaks == 21) %>%
left_join(meth, by = c("year" = "Year")) %>%
select(year, Probability, Length) %>%
mutate(days = Probability*Length)
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogaW5saW5lCi0tLQoKRGF0YSBmcm9tIDE5NjAgdG8gMjAxMCByYWluaW5nIGRheXMsIHRoZXJlIGFyZSA0MTY4IG9ic2VydmF0aW9ucyAxMyBkYXlzIG9mIHdoaWNoIGFyZSByYWluIGV4Y2VlZHMgMTAwMCAoMTBjbSkKCkNoYW5nZSBvZiBleGNlZWRzIDEwMDAgaXMgMC4wMDMxMTkwMDIgMC4zMSUgY2hhbmdlCgpgYGB7cn0KCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KHJlYWRyKQpsaWJyYXJ5KGRyYykKbGlicmFyeShyZXNoYXBlMikKCmluc3RhbGwucGFja2FnZXMoIm1vbWVudHMiKSAgICAgICAgICAgICAgICAKbGlicmFyeSgibW9tZW50cyIpICMgZm9yIHNrZXduZXNzIG9mIGRhdGEKYGBgCgojIERhdGEKCmBgYHtyfQpyZWFkLmNzdigiLi9yYXcvcHJjcC5jc3YiKSAlPiUKICBtdXRhdGUoZGF0ZSA9IGFzLkRhdGUoZGF0ZSkpICU+JQogIG11dGF0ZSh5ZWFyID0gZm9ybWF0KGRhdGUsICIlWSIpKSAlPiUKICBmaWx0ZXIocHJjcCAhPSAwKSAlPiUKICBkcGx5cjo6c2VsZWN0KGRhdGUsIHllYXIsIHByY3ApIC0+IGRhdGEKCmRhdGEgJT4lCiAgZmlsdGVyKHllYXIgPD0gMjAyMCkgJT4lCiAgcHVsbChwcmNwKSAlPiUKICBoaXN0KGJyZWFrcyA9IDEwMDAsIHBsb3QgPSBGQUxTRSkgLT4gaGlzdApgYGAKCiMgTW9kbGluZyBEaXN0cmlidXRpb24KClVzZSBkYXRhIGZyb20gMTk2MCAtIDIwMTAgZm9yIG1vZGVsaW5nIHByZWNpcGl0YXRpb24KCmBgYHtyfQpoaXN0JGNvdW50cyAlPiUgc3VtKCkgLT4gc3VtCnggPC0gaGlzdCRtaWRzIAp5IDwtIGhpc3QkY291bnRzL3N1bSAKCmRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKYGBgCgpgYGB7cn0KbW9kZWwgPC0gZHJtKGRhdGEgPSBkZW5zLCB5IH4gbG9nKHgpLCBmY3QgPSBFWEQuMigpKQpzdW1tYXJ5KG1vZGVsKQpkID0gMC4xNTU2NzQ5OQplID0gMS4yNzQyNDgzMgoKayA8LSAxL2UvZApgYGAKCmBgYHtyfQpwbG90KHggPSBkZW5zJHgsIHkgPSBmaXR0ZWQobW9kZWwpLCB0eXBlID0gImwiLCBjb2wgPSAiYmx1ZSIsIHhsaW0gPSBjKDAsIDMwMCkpCmxpbmVzKHggPSBkZW5zJHgsIHkgPSBkZW5zJHksIHR5cGUgPSAicCIsIGNvbCA9IGFscGhhKCJibGFjayIsIDAuNSkpCmBgYAoKYGBge3J9Cm1vZC5leGQgPC0gZHJtKHkgfiB4LCBmY3QgPSBFWEQuMigpLCBkYXRhID0gZGVucykKYGBgCgpgYGB7cn0KcGxvdCh4ID0gZGVucyR4LCB5ID0gZml0dGVkKG1vZC5leGQpLCB0eXBlID0gImwiLCBjb2wgPSAiYmx1ZSIsIHhsaW0gPSBjKDAsMzAwKSkKbGluZXMoeCA9IGRlbnMkeCwgeSA9IGRlbnMkeSwgdHlwZSA9ICJwIiwgY29sID0gYWxwaGEoImJsYWNrIiwgMC41KSkKYGBgCgpgYGB7cn0Kc3VtbWFyeShtb2QuZXhkKQoKYGBgCgpgYGB7cn0KMSAtIHBleHAobG9nKDEwMDApKmssIGQpCgoxL2QKCmRhdGEgJT4lCiAgZmlsdGVyKHllYXIgPD0gMjAxMCkgJT4lCiAgcHVsbChwcmNwKSAlPiUKICBtZWFuKCkgJT4lCiAgbG9nKCkqawpgYGAKCjAuMzElIGNoYW5nZSB2ZXJzdXMgMC40NDIlCgojIExhYnM6CgojIyBBIFNpbmdsZSBZZWFycwoKYGBge3J9CnlyID0gMTk2MApkYXRhICU+JQogIGZpbHRlcigoeWVhciA+PSB5cikgJiAoeWVhciA8PSB5ciArIDkpKSAlPiUKICBwdWxsKHByY3ApICU+JQogIGhpc3QoYnJlYWtzID0gMTAwMCwgcGxvdCA9IEZBTFNFKSAtPiBoaXN0Cmhpc3QkY291bnRzICU+JSBzdW0oKSAtPiBzdW0KeCA8LSBoaXN0JG1pZHMgCnkgPC0gaGlzdCRjb3VudHMvc3VtIApkYXRhLmZyYW1lKHgsIHkpIC0+IGRlbnMgCm1vZCA8LSBkcm0oZGF0YSA9IGRlbnMsIHkgfiBsb2coeCksIGZjdCA9IEVYRC4yKCkpCmQgPC0gbW9kJHBhcm1NYXRbMV0KZSA8LSBtb2QkcGFybU1hdFsyXQprIDwtIDEvZS9kCjEgLSBwZXhwKGxvZygxMDAwKSprLCBkKQpgYGAKCmBgYHtyfQp5ciA9IDE5NzAKZGF0YSAlPiUKICBmaWx0ZXIoKHllYXIgPj0geXIpICYgKHllYXIgPD0geXIgKyA5KSkgJT4lCiAgcHVsbChwcmNwKSAlPiUKICBoaXN0KGJyZWFrcyA9IDEwMDAsIHBsb3QgPSBGQUxTRSkgLT4gaGlzdApoaXN0JGNvdW50cyAlPiUgc3VtKCkgLT4gc3VtCnggPC0gaGlzdCRtaWRzIAp5IDwtIGhpc3QkY291bnRzL3N1bSAKZGF0YS5mcmFtZSh4LCB5KSAtPiBkZW5zIAptb2QgPC0gZHJtKGRhdGEgPSBkZW5zLCB5IH4gbG9nKHgpLCBmY3QgPSBFWEQuMigpKQpkIDwtIG1vZCRwYXJtTWF0WzFdCmUgPC0gbW9kJHBhcm1NYXRbMl0KayA8LSAxL2UvZAoxIC0gcGV4cChsb2coMTAwMCkqaywgZCkKYGBgCgpgYGB7cn0KeXIgPSAxOTgwCmRhdGEgJT4lCiAgZmlsdGVyKCh5ZWFyID49IHlyKSAmICh5ZWFyIDw9IHlyICsgOSkpICU+JQogIHB1bGwocHJjcCkgJT4lCiAgaGlzdChicmVha3MgPSAxMDAwLCBwbG90ID0gRkFMU0UpIC0+IGhpc3QKaGlzdCRjb3VudHMgJT4lIHN1bSgpIC0+IHN1bQp4IDwtIGhpc3QkbWlkcyAKeSA8LSBoaXN0JGRlbnNpdHkgCmRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKbW9kIDwtIGRybShkYXRhID0gZGVucywgeSB+IGxvZyh4KSwgZmN0ID0gRVhELjIoKSkKZCA8LSBtb2QkcGFybU1hdFsxXQplIDwtIG1vZCRwYXJtTWF0WzJdCmsgPC0gMS9lL2QKMSAtIHBleHAobG9nKDEwMDApKmssIGQpCmBgYAoKIyMgUmlzayBUcmVuZHMKCiMjIFJpc2tzIFRocm91Z2ggb3V0IFllYXJzCgpgYGB7cn0KbXRyIDwtIGRhdGEuZnJhbWUoeWVhciA9IFlFQVIpCmxpc3QgPC0gYygpCgpmb3IgKGkgaW4gWUVBUil7CiAgICBkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0gaSkgJiAoeWVhciA8PSBpICsgOSkpICU+JQogICAgICBwdWxsKHByY3ApICU+JQogICAgICBoaXN0KGJyZWFrcyA9IDQwMCwgcGxvdCA9IEZBTFNFKSAtPiBoaXN0CiAgICBoaXN0JGNvdW50cyAlPiUgc3VtKCkgLT4gc3VtCiAgICB4IDwtIGhpc3QkbWlkcyAKICAgIHkgPC0gaGlzdCRjb3VudHMvc3VtIAogICAgZGF0YS5mcmFtZSh4LCB5KSAtPiBkZW5zIAogICAgbW9kIDwtIGRybShkYXRhID0gZGVucywgeSB+IGxvZyh4KSwgZmN0ID0gRVhELjIoKSkKICAgIGQgPC0gbW9kJHBhcm1NYXRbMV0KICAgIGUgPC0gbW9kJHBhcm1NYXRbMl0KICAgIGsgPC0gMS9lL2QKICAgIHAgPC0gMSAtIHBleHAobG9nKDEwMDApKmssIGQpCiAgICBsaXN0IDwtIGMobGlzdCwgcCkKfQpwbG90KDE6bGVuZ3RoKGxpc3QpLCBsaXN0KQpgYGAKCmBgYHtyfQoKYGBgCgojIyBVc2Ugc21hbGxlciBicmVha3MKCkZpbmQgdGhlIHJpZ2h0IGJyZWFrIG51bWJlcnMKCmBgYHtyfQpZRUFSIDwtIGMoMTk2MCwgMTk4MCwgMjAwMCkKQlJFQUtTIDwtIGMoMTksIDIwLCAyMSwgMjIpCmBgYAoKYGBge3J9Cmxlbmxpc3QgPC0gYygpCmt1bGlzdCA8LSBjKCkKZm9yKHlyIGluIFlFQVIpewpkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0geXIpICYgKHllYXIgPD0geXIgKyAxOSkpICU+JQogICAgICBkcGx5cjo6c2VsZWN0KHllYXIsIHByY3ApICU+JSAKICAgICAgcHVsbChwcmNwKSAlPiUKICAgICAgbGVuZ3RoKCkgLT4gbgogIGxlbmxpc3QgPC0gYyhsZW5saXN0LCBuKQpkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0geXIpICYgKHllYXIgPD0geXIgKyAxOSkpICU+JQogICAgICBkcGx5cjo6c2VsZWN0KHllYXIsIHByY3ApICU+JSAKICAgICAgcHVsbChwcmNwKSAlPiUKICAgICAga3VydG9zaXMoKSAtPiBrdQogIGt1bGlzdCA8LSBjKGt1bGlzdCwga3UpCn0KU3F1YXJlX3Jvb3QgPSBjKCkKU3R1cmdlcyA9IGMoKQpSaWNlID0gYygpCkRvbm5hID0gYygpCmZvcihuIGluIGxlbmxpc3QpewogIHNxID0gY2VpbGluZyhuXjAuNSkKICBzZyA9IGNlaWxpbmcobG9nKG4sIDIpKSArIDEKICByYyA9IGNlaWxpbmcoMipuXigxLzMpKQogIGlnaCA9ICg2KihuIC0gMikvKG4gKzEpLyhuICszKSleMC41CiAgZG4gPSBsb2coYWJzKDEgKyBrdWxpc3Rbd2hpY2gobGVubGlzdCAlaW4lIG4pXS9pZ2gpICwgMikgKyBsb2cobiwgMikgKyAxClNxdWFyZV9yb290ID0gYyhTcXVhcmVfcm9vdCwgc3EpClN0dXJnZXMgPSBjKFN0dXJnZXMsIHNnKQpSaWNlID0gYyhSaWNlLCByYykKRG9ubmEgPSBjKERvbm5hLCBkbikKfQpkYXRhLmZyYW1lKCBZZWFyID0gWUVBUiwKICAgICAgICAgIExlbmd0aCA9IGxlbmxpc3QsCiAgICAgICAgICAgU3F1YXJlX3Jvb3QgPSBTcXVhcmVfcm9vdCwgCiAgICAgICAgICAgU3R1cmdlcyA9IFN0dXJnZXMsCiAgICAgICAgICAgUmljZSA9IFJpY2UsIAogICAgICAgICAgIERvbm5hID0gRG9ubmEKICAgICAgICAgICApIC0+IG1ldGgKYGBgCgpgYGB7cn0KbWV0aApgYGAKCiMgV29yayBGbG93CgpJbiBhIDEwIHllYXJzIGludGVydmFsIGNvbXB1dGUgcmlza3Mgb2YgcHJlY2lwaXRhdGlvbiBoaWdoZXIgdGhhbiAxMGNtIChwcmNwIFw+IDEwMDApIGJ5IHVzaW5nIGRpZmZlcmVudCBicmVha3Mgb2YgaGlzdG9ncmFtIGZyb20gdXNnaW4gMjUwIHRvIDEwMDAKCkNob29zaW5nIGJpbiBoYXZlIG1hbnkgY2hvaWNlcyBhY3R1YWxseQoKYGBge3J9CllFQVIgPC0gYygxOTYwLCAxOTgwLCAyMDAwKQpCUkVBS1MgPC0gYygxOSwgMjAsIDIxLCAyMiwgMzAsIDQwLCA1MCwgNjAsIDcwLCA4MCwgOTApCmBgYAoKYGBge3J9CmZvcihuIGluIEJSRUFLUykKICBkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0gMTk2MCkgJiAoeWVhciA8PSAxOTY5KSkgJT4lCiAgICAgIHB1bGwocHJjcCkgJT4lCiAgICAgIGhpc3QoYnJlYWtzID0gbikKYGBgCgpgYGB7cn0Ka3VydG9zaXMoZGF0YSRwcmNwKQpgYGAKCmBgYHtyfQpsaXN0IDwtIGMoKQpsaXN0QiA8LSBjKCkKWWVhciA8LSBjKCkKZGxpc3QgPC0gYygpCmVsaXN0IDwtIGMoKQprbGlzdCA8LSBjKCkKZm9yKG4gaW4gQlJFQUtTKXsKICBmb3IoaSBpbiBZRUFSKXsKICAgIGRhdGEgJT4lCiAgICAgIGZpbHRlcigoeWVhciA+PSBpKSAmICh5ZWFyIDw9IGkgKyAxOSkpICU+JQogICAgICBwdWxsKHByY3ApICU+JQogICAgICBoaXN0KGJyZWFrcyA9IG4sIHBsb3QgPSBGQUxTRSkgLT4gaGlzdAogICAgaGlzdCRjb3VudHMgJT4lIHN1bSgpIC0+IHN1bQogICAgeCA8LSBoaXN0JG1pZHMgCiAgICB5IDwtIGhpc3QkZGVuc2l0eSAKICAgIGRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKICAgIG1vZCA8LSBkcm0oZGF0YSA9IGRlbnMsIHkgfiBsb2coeCksIGZjdCA9IEVYRC4yKCkpCiAgICBkIDwtIG1vZCRwYXJtTWF0WzFdCiAgICBlIDwtIG1vZCRwYXJtTWF0WzJdCiAgICBrIDwtIDEvZS9kCiAgICBwIDwtIDEgLSBwZXhwKGxvZygxNTAwKSprLCBkKQogICAgbGlzdCA8LSBjKGxpc3QsIHApCiAgICBsaXN0QiA8LSBjKGxpc3RCLCBuKQogICAgWWVhciA8LSBjKFllYXIsIGkpCiAgICBkbGlzdCA8LSBjKGRsaXN0LCBkKQogICAgZWxpc3QgPC0gYyhlbGlzdCwgZSkKICAgIGtsaXN0IDwtIGMoa2xpc3QsIGspCiAgfX0KbXRyIDwtIGRhdGEuZnJhbWUoeWVhciA9IFllYXIsCiAgICAgICAgICAgICAgICAgICAgQnJlYWtzID0gbGlzdEIsIAogICAgICAgICAgICAgICAgICAgICAgUHJvYmFiaWxpdHkgPSBsaXN0LCAKICAgICAgICAgICAgICAgICAgcGFyYW1ldGVyX2QgPSBkbGlzdCwKICAgICAgICAgICAgICAgICAgcGFyYW1ldGVyX2UgPSBlbGlzdCwKICAgICAgICAgICAgICAgICAgcGFyYW1ldGVyX2sgPSBrbGlzdCkKYGBgCgojIFZpc3VsaXNhdGlvbgoKYGBge3J9CmxpYnJhcnkod2VzYW5kZXJzb24pCm10ciAlPiUKICBtdXRhdGUoQnJlYWtzID0gYXMuZmFjdG9yKEJyZWFrcykpICU+JQogIGdncGxvdChhZXMoeCA9IHllYXIsIHkgPSBQcm9iYWJpbGl0eSwgY29sb3IgPSBCcmVha3MpKSArIAogICAgZ2VvbV9wb2ludCgpICsgCiAgICBnZW9tX2xpbmUoKSArIAogICAgc2NhbGVfY29sb3JfYnJld2VyKHBhbGV0dGUgPSAiQmx1ZXMiKSArCiAgdGhlbWUocGFuZWwuYm9yZGVyID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5tYWpvciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWlub3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KGZpbGwgPSAiI2ZhZmFmYSIpLAogICAgICAgICAgcGxvdC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KGZpbGwgPSAid2hpdGUiKQogICAgICAgICkgKyAKICAgIHlsYWIoIlJpc2tzIG9mIFJhaW4gRXhjZWVkcyAxMGNtIikgKyAKICAgIGdndGl0bGUoIlJpc2tzIG9mIFJpYW4gRXhjZWVkaW5nIDE1IGNtIGluY3JlYXNlIGV2ZXkgMiBkZWNhZGVzIikgKwogICAgeGxhYigiVHdvIERlY2FkZXMgVGltZSBGcm9tIF8gIikgLT4gZwojbGlicmFyeShwbG90bHkpCmdncGxvdGx5KGcpCmBgYAoKYGBge3J9Cm10ciAlPiUKICBtdXRhdGUoeWVhciA9IGFzLmZhY3Rvcih5ZWFyKSkgJT4lCiAgbXV0YXRlKEJyZWFrcyA9IGFzLmZhY3RvcihCcmVha3MpKSAlPiUKICBnZ3Bsb3QoYWVzKHkgPSBwYXJhbWV0ZXJfZCwgeCA9IHBhcmFtZXRlcl9lLCBjb2xvciA9IHllYXIpKSArIAogICAgZ2VvbV9qaXR0ZXIoYWxwaGEgPSAwLjgpICsgCiAgICBzY2FsZV9jb2xvcl9icmV3ZXIocGFsZXR0ZSA9ICJSZWRzIikgKwogIHRoZW1lKHBhbmVsLmJvcmRlciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1pbm9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIiNmYWZhZmEiKSwKICAgICAgICAgIHBsb3QuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIndoaXRlIikKICAgICAgICApICsgCiAgICB4bGFiKCJwYXJhbWV0ZXIgZTogc2tld25lc3Mgb2YgdGhlIGN1cnZlIikgKyAKICAgIHlsYWIoInBhcmFtZXRlciBkOiBoZWlnaHQgb2YgdGhlIGN1cnZlIikgCiAgICAKYGBgCgpQYXJhbWV0ZXIgZCBpbmZsdWVuY2UgdGhlIGhlaWdodCBvZiBwZGYsIHBhcmFtZXRlciBlIGluZmx1ZW5jZSB0aGUgbGVuZ2h0aCBvZiBwZGYuLi4gdGhpcyBzdWdnZXN0IGluIGdlbmVyYWwsIHByZWNpcGl0YXRpb24gZGlzdHJpYnV0aW9uIGN1cmUgYXJlIGdldHRpbmcgZmxhdHRlciBhbmQKCmBgYHtyfQpwbG90X2x5KHg9IG10ciRwYXJhbWV0ZXJfZSwgeT0gbXRyJHBhcmFtZXRlcl9kLCB6PW10ciRCcmVha3MsIHR5cGU9InNjYXR0ZXIzZCIsIG1vZGU9Im1hcmtlcnMiLCBjb2xvcj1tdHIkeWVhcikKYGBgCgpgYGB7cn0KbXRyICU+JQogIGRwbHlyOjpzZWxlY3QoeWVhciwgQnJlYWtzLCBwYXJhbWV0ZXJfZCwgcGFyYW1ldGVyX2UpICU+JQogIG11dGF0ZShCcmVha3MgPSBhcy5mYWN0b3IoQnJlYWtzKSkgJT4lCiAgbWVsdChpZCA9IGMoInllYXIiLCAiQnJlYWtzIiksIHZhcmlhYmxlLm5hbWUgPSAicGFyYW1ldGVyIikgJT4lCiAgZ2dwbG90KGFlcyh4ID0geWVhciwgeSA9IHZhbHVlKSkgKyAKICAgIGdlb21fbGluZShhZXMoeCA9IHllYXIsIHkgPSB2YWx1ZSwgbGluZXR5cGUgPSBwYXJhbWV0ZXIsIGNvbG9yID0gQnJlYWtzKSkgKwogICAgc2NhbGVfY29sb3JfYnJld2VyKHBhbGV0dGUgPSAiUHVycGxlcyIpICsKICB0aGVtZShwYW5lbC5ib3JkZXIgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1ham9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICIjZmFmYWZhIiksCiAgICAgICAgICBwbG90LmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICJ3aGl0ZSIpCiAgICAgICAgKQpgYGAKCmBgYHtyfQptdHIgJT4lCiAgbXV0YXRlKEJyZWFrcyA9IGFzLmZhY3RvcihCcmVha3MpKSAlPiUKICBnZ3Bsb3QoYWVzKHggPSB5ZWFyLCB5ID0gcGFyYW1ldGVyX2ssIGNvbG9yID0gQnJlYWtzKSkgKyAKICAgIGdlb21fcG9pbnQoKSArIAogICAgZ2VvbV9saW5lKCkgKyAKICAgIHNjYWxlX2NvbG9yX2JyZXdlcihwYWxldHRlID0gIkJsdWVzIikgKwogIHRoZW1lKHBhbmVsLmJvcmRlciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1pbm9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIiNmYWZhZmEiKSwKICAgICAgICAgIHBsb3QuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIndoaXRlIikKICAgICAgICApCmBgYAoKY2FsY3VsYXRlIGRpZmZlcmVudAoKYGBge3J9Cmxpc3QgPC0gYygpCmxpc3RCIDwtIGMoKQpZZWFyIDwtIGMoKQpkbGlzdCA8LSBjKCkKZWxpc3QgPC0gYygpCmZvcihuIGluIEJSRUFLUyl7CiAgZm9yKGkgaW4gWUVBUil7CiAgICBkYXRhICU+JQogICAgICBmaWx0ZXIoKHllYXIgPj0gaSkgJiAoeWVhciA8PSBpICsgOSkpICU+JQogICAgICBwdWxsKHByY3ApICU+JQogICAgICBoaXN0KGJyZWFrcyA9IG4sIHBsb3QgPSBGQUxTRSkgLT4gaGlzdAogICAgaGlzdCRjb3VudHMgJT4lIHN1bSgpIC0+IHN1bQogICAgeCA8LSBoaXN0JG1pZHMgCiAgICB5IDwtIGhpc3QkY291bnRzL3N1bSAKICAgIGRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKICAgIG1vZCA8LSBkcm0oZGF0YSA9IGRlbnMsIHkgfiBsb2coeCksIGZjdCA9IEVYRC4yKCkpCiAgICBkIDwtIG1vZCRwYXJtTWF0WzFdCiAgICBlIDwtIG1vZCRwYXJtTWF0WzJdCiAgICBrIDwtIDEvZS9kCiAgICBwIDwtIDEgLSBwZXhwKGxvZyg2MDApKmssIGQpCiAgICBsaXN0IDwtIGMobGlzdCwgcCkKICAgIGxpc3RCIDwtIGMobGlzdEIsIG4pCiAgICBZZWFyIDwtIGMoWWVhciwgaSkKICAgIGRsaXN0IDwtIGMoZGxpc3QsIGQpCiAgICBlbGlzdCA8LSBjKGVsaXN0LCBlKQogIH19Cm10ciA8LSBkYXRhLmZyYW1lKHllYXIgPSBZZWFyLAogICAgICAgICAgICAgICAgICAgIEJyZWFrcyA9IGxpc3RCLCAKICAgICAgICAgICAgICAgICAgICAgIFByb2JhYmlsaXR5ID0gbGlzdCwgCiAgICAgICAgICAgICAgICAgIHBhcmFtZXRlcl9kID0gZGxpc3QsCiAgICAgICAgICAgICAgICAgIHBhcmFtZXRlcl9lID0gZWxpc3QpCgpgYGAKCmBgYHtyfQptdHIgJT4lCiAgbXV0YXRlKEJyZWFrcyA9IGFzLmZhY3RvcihCcmVha3MpKSAlPiUKICBnZ3Bsb3QoYWVzKHggPSB5ZWFyLCB5ID0gUHJvYmFiaWxpdHksIGNvbG9yID0gQnJlYWtzKSkgKyAKICAgIGdlb21fcG9pbnQoKSArIAogICAgZ2VvbV9saW5lKCkgKyAKICAgIHNjYWxlX2NvbG9yX2JyZXdlcihwYWxldHRlID0gIkJsdWVzIikgKwogIHRoZW1lKHBhbmVsLmJvcmRlciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1pbm9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIiNmYWZhZmEiKSwKICAgICAgICAgIHBsb3QuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdChmaWxsID0gIndoaXRlIikKICAgICAgICApCmBgYAoKYGBge3J9CmcgPC0gZGYgJT4lCiAgZ2dwbG90KGFlcyh4ID0gZGF0ZSwgeSA9IHByY3ApKSArIAogIGdlb21fbGluZShjb2xvciA9ICJibHVlIikgKyAKICB0aGVtZShwYW5lbC5ib3JkZXIgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1ham9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICIjZmFmYWZhIiksCiAgICAgICAgICBwbG90LmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICJ3aGl0ZSIpCiAgICAgICAgKSArIAogIGdlb21fbGluZShhZXMoeCA9IGRhdGUsIHkgPSByZXAoMTUwMCwgbGVuZ3RoKGRhdGUpKSksIGNvbG9yID0gImxpZ2h0Ymx1ZSIpICsgCiAgeWxhYigiUHJlY2lwaXRhdGlvbiAoaW4gdGVudGggb2YgbW0iKQpnZ3Bsb3RseShnKQpgYGAKCmBgYHtyfQpkYXRhICU+JQogIGZpbHRlcihwcmNwID49IDEwMDApCmBgYAoKYGBge3J9Cm10ciAlPiUKICBmaWx0ZXIoQnJlYWtzID09IDIxKSAlPiUKICBsZWZ0X2pvaW4obWV0aCwgYnkgPSBjKCJ5ZWFyIiA9ICJZZWFyIikpICU+JQogIHNlbGVjdCh5ZWFyLCBQcm9iYWJpbGl0eSwgTGVuZ3RoKSAlPiUKICBtdXRhdGUoZGF5cyA9IFByb2JhYmlsaXR5Kkxlbmd0aCkKYGBgCg==